home *** CD-ROM | disk | FTP | other *** search
/ AP Professional Graphics CD-ROM Library / AP Professional Graphics CD-ROM Library.iso / pc / unix / appendix / gemsi / root3n4.c < prev    next >
Encoding:
C/C++ Source or Header  |  1992-09-15  |  4.7 KB  |  245 lines
  1. /*
  2.  *  Roots3And4.c
  3.  *
  4.  *  Utility functions to find cubic and quartic roots,
  5.  *  coefficients are passed like this:
  6.  *
  7.  *      c[0] + c[1]*x + c[2]*x^2 + c[3]*x^3 + c[4]*x^4 = 0
  8.  *
  9.  *  The functions return the number of non-complex roots and
  10.  *  put the values into the s array.
  11.  *
  12.  *  Author:         Jochen Schwarze (schwarze@isa.de)
  13.  *
  14.  *  Jan 26, 1990    Version for Graphics Gems
  15.  *  Oct 11, 1990    Fixed sign problem for negative q's in SolveQuartic
  16.  *                  (reported by Mark Podlipec),
  17.  *                  Old-style function definitions,
  18.  *                  IsZero() as a macro
  19.  *  Nov 23, 1990    Some systems do not declare acos() and cbrt() in
  20.  *                  <math.h>, though the functions exist in the library.
  21.  *                  If large coefficients are used, EQN_EPS should be
  22.  *                  reduced considerably (e.g. to 1E-30), results will be
  23.  *                  correct but multiple roots might be reported more
  24.  *                  than once.
  25.  */
  26.  
  27. #include    <math.h>
  28. #ifndef M_PI
  29. #define M_PI          3.14159265358979323846
  30. #endif
  31. extern double   sqrt(), cbrt(), cos(), acos();
  32.  
  33. /* epsilon surrounding for near zero values */
  34.  
  35. #define     EQN_EPS     1e-9
  36. #define        IsZero(x)    ((x) > -EQN_EPS && (x) < EQN_EPS)
  37.  
  38. #ifdef NOCBRT
  39. #define     cbrt(x)     ((x) > 0.0 ? pow((double)(x), 1.0/3.0) : \
  40.                           ((x) < 0.0 ? -pow((double)-(x), 1.0/3.0) : 0.0))
  41. #endif
  42.  
  43. int SolveQuadric(c, s)
  44.     double c[ 3 ];
  45.     double s[ 2 ];
  46. {
  47.     double p, q, D;
  48.  
  49.     /* normal form: x^2 + px + q = 0 */
  50.  
  51.     p = c[ 1 ] / (2 * c[ 2 ]);
  52.     q = c[ 0 ] / c[ 2 ];
  53.  
  54.     D = p * p - q;
  55.  
  56.     if (IsZero(D))
  57.     {
  58.     s[ 0 ] = - p;
  59.     return 1;
  60.     }
  61.     else if (D < 0)
  62.     {
  63.     return 0;
  64.     }
  65.     else if (D > 0)
  66.     {
  67.     double sqrt_D = sqrt(D);
  68.  
  69.     s[ 0 ] =   sqrt_D - p;
  70.     s[ 1 ] = - sqrt_D - p;
  71.     return 2;
  72.     }
  73. }
  74.  
  75.  
  76. int SolveCubic(c, s)
  77.     double c[ 4 ];
  78.     double s[ 3 ];
  79. {
  80.     int     i, num;
  81.     double  sub;
  82.     double  A, B, C;
  83.     double  sq_A, p, q;
  84.     double  cb_p, D;
  85.  
  86.     /* normal form: x^3 + Ax^2 + Bx + C = 0 */
  87.  
  88.     A = c[ 2 ] / c[ 3 ];
  89.     B = c[ 1 ] / c[ 3 ];
  90.     C = c[ 0 ] / c[ 3 ];
  91.  
  92.     /*  substitute x = y - A/3 to eliminate quadric term:
  93.     x^3 +px + q = 0 */
  94.  
  95.     sq_A = A * A;
  96.     p = 1.0/3 * (- 1.0/3 * sq_A + B);
  97.     q = 1.0/2 * (2.0/27 * A * sq_A - 1.0/3 * A * B + C);
  98.  
  99.     /* use Cardano's formula */
  100.  
  101.     cb_p = p * p * p;
  102.     D = q * q + cb_p;
  103.  
  104.     if (IsZero(D))
  105.     {
  106.     if (IsZero(q)) /* one triple solution */
  107.     {
  108.         s[ 0 ] = 0;
  109.         num = 1;
  110.     }
  111.     else /* one single and one double solution */
  112.     {
  113.         double u = cbrt(-q);
  114.         s[ 0 ] = 2 * u;
  115.         s[ 1 ] = - u;
  116.         num = 2;
  117.     }
  118.     }
  119.     else if (D < 0) /* Casus irreducibilis: three real solutions */
  120.     {
  121.     double phi = 1.0/3 * acos(-q / sqrt(-cb_p));
  122.     double t = 2 * sqrt(-p);
  123.  
  124.     s[ 0 ] =   t * cos(phi);
  125.     s[ 1 ] = - t * cos(phi + M_PI / 3);
  126.     s[ 2 ] = - t * cos(phi - M_PI / 3);
  127.     num = 3;
  128.     }
  129.     else /* one real solution */
  130.     {
  131.     double sqrt_D = sqrt(D);
  132.     double u = cbrt(sqrt_D - q);
  133.     double v = - cbrt(sqrt_D + q);
  134.  
  135.     s[ 0 ] = u + v;
  136.     num = 1;
  137.     }
  138.  
  139.     /* resubstitute */
  140.  
  141.     sub = 1.0/3 * A;
  142.  
  143.     for (i = 0; i < num; ++i)
  144.     s[ i ] -= sub;
  145.  
  146.     return num;
  147. }
  148.  
  149.  
  150. int SolveQuartic(c, s)
  151.     double c[ 5 ]; 
  152.     double s[ 4 ];
  153. {
  154.     double  coeffs[ 4 ];
  155.     double  z, u, v, sub;
  156.     double  A, B, C, D;
  157.     double  sq_A, p, q, r;
  158.     int     i, num;
  159.  
  160.     /* normal form: x^4 + Ax^3 + Bx^2 + Cx + D = 0 */
  161.  
  162.     A = c[ 3 ] / c[ 4 ];
  163.     B = c[ 2 ] / c[ 4 ];
  164.     C = c[ 1 ] / c[ 4 ];
  165.     D = c[ 0 ] / c[ 4 ];
  166.  
  167.     /*  substitute x = y - A/4 to eliminate cubic term:
  168.     x^4 + px^2 + qx + r = 0 */
  169.  
  170.     sq_A = A * A;
  171.     p = - 3.0/8 * sq_A + B;
  172.     q = 1.0/8 * sq_A * A - 1.0/2 * A * B + C;
  173.     r = - 3.0/256*sq_A*sq_A + 1.0/16*sq_A*B - 1.0/4*A*C + D;
  174.  
  175.     if (IsZero(r))
  176.     {
  177.     /* no absolute term: y(y^3 + py + q) = 0 */
  178.  
  179.     coeffs[ 0 ] = q;
  180.     coeffs[ 1 ] = p;
  181.     coeffs[ 2 ] = 0;
  182.     coeffs[ 3 ] = 1;
  183.  
  184.     num = SolveCubic(coeffs, s);
  185.  
  186.     s[ num++ ] = 0;
  187.     }
  188.     else
  189.     {
  190.     /* solve the resolvent cubic ... */
  191.  
  192.     coeffs[ 0 ] = 1.0/2 * r * p - 1.0/8 * q * q;
  193.     coeffs[ 1 ] = - r;
  194.     coeffs[ 2 ] = - 1.0/2 * p;
  195.     coeffs[ 3 ] = 1;
  196.  
  197.     (void) SolveCubic(coeffs, s);
  198.  
  199.     /* ... and take the one real solution ... */
  200.  
  201.     z = s[ 0 ];
  202.  
  203.     /* ... to build two quadric equations */
  204.  
  205.     u = z * z - r;
  206.     v = 2 * z - p;
  207.  
  208.     if (IsZero(u))
  209.         u = 0;
  210.     else if (u > 0)
  211.         u = sqrt(u);
  212.     else
  213.         return 0;
  214.  
  215.     if (IsZero(v))
  216.         v = 0;
  217.     else if (v > 0)
  218.         v = sqrt(v);
  219.     else
  220.         return 0;
  221.  
  222.     coeffs[ 0 ] = z - u;
  223.     coeffs[ 1 ] = q < 0 ? -v : v;
  224.     coeffs[ 2 ] = 1;
  225.  
  226.     num = SolveQuadric(coeffs, s);
  227.  
  228.     coeffs[ 0 ]= z + u;
  229.     coeffs[ 1 ] = q < 0 ? v : -v;
  230.     coeffs[ 2 ] = 1;
  231.  
  232.     num += SolveQuadric(coeffs, s + num);
  233.     }
  234.  
  235.     /* resubstitute */
  236.  
  237.     sub = 1.0/4 * A;
  238.  
  239.     for (i = 0; i < num; ++i)
  240.     s[ i ] -= sub;
  241.  
  242.     return num;
  243. }
  244.  
  245.